#include "StVenKir.hpp"

// variational formulation language
#include <life/lifevf/vf.hpp>

namespace Life
{


  template<int Dim, int Order, template<uint16_type,uint16_type,uint16_type> class Entity>
  void
  StVenantKirchhoff<Dim, Order, Entity>::updateJacobian( const vector_ptrtype& X, sparse_matrix_ptrtype& J)
  {

    using namespace Life::vf;

    mesh_ptrtype mesh = M_Xh->mesh();

    element_struct_type u( M_Xh, "u" ); u = *X;
    element_struct_type v( M_Xh, "v" );



    //--------------------------------------------------------------------------------------------------//

    //coeff ...
    //   double lambda = this->vm()["rho"].template as<double>();
    //  double nu = this->vm()["rho"].template as<double>();


    //Matrice Identitée
    AUTO( Id , oneX()*trans(oneX()) + oneY()*trans(oneY()) );

    AUTO( Fv , Id + gradv(u) );
    AUTO( F , Id + grad(u) );
    AUTO( dF, gradt(u));
    AUTO( dFF, grad(u));

    AUTO( Ev, 0.5*(gradv(u)+trans(gradv(u)) ) + 0.5*trans(gradv(u))*gradv(u) );
    AUTO( E, 0.5*(grad(u)+trans(grad(u)) ) + 0.5*trans(grad(u))*grad(u) );
    AUTO( dE, 0.5*(trans(Fv)*dF + trans(dF)*Fv) );
    AUTO( dEE, 0.5*(trans(Fv)*dFF + trans(dFF)*Fv) );

    AUTO( Sv, lambda*trace(Ev)*Id + 2*mu*Ev );
    AUTO( S, lambda*trace(E)*Id + 2*mu*E );
    AUTO( dS, lambda*trace(dE)*Id + 2*mu*dE );
    AUTO( dSS, lambda*trace(dEE)*Id + 2*mu*dEE );

    //--------------------------------------------------------------------------------------------------//
    form2( M_Xh, M_Xh, J, _init=true )  =
      integrate (elements(mesh), _Q<4*Order>(),
		 trace( (dF*Sv + Fv*dS)*trans(grad(v)) ) );

    //Attention les termes de symetrisation ne sont pas necessaire mais accelere la convergence
    // Prendre le parametre de stabilisation d'environ 1e6!!!
    double penalisation_bc=this->application()->vm()["bccoeff"].template as<double>();
#if 1
    form2( M_Xh, M_Xh, J )  +=
      integrate( markedfaces(mesh,mesh->markerName("Fixe")), _Q<4*Order>(),
                 //terme classique
                 - trans((dF*Sv + Fv*dS)*N())*id(v)
                 //terme de symetrisation
                 - trans((dFF*Sv + Fv*dSS)*N())*idt(u)
                 //            - trans((dF*Sv + Fv*dS)*N())*idt(u)
                 //            - trans((F*S)*N())*idt(u)
                 + penalisation_bc*trans(idt(u))*id(v)/hFace() );
#endif
#if 1
    form2( M_Xh, M_Xh, J )  +=
      integrate( markedfaces(mesh,mesh->markerName("Force")), _Q<4*Order>(),
                 //terme classique
                 - trans((dF*Sv + Fv*dS)*N())*id(v)
                 //terme de symetrisation
                 - trans((dFF*Sv + Fv*dSS)*N())*idt(u)
                 //            - trans((dF*Sv + Fv*dS)*N())*idt(u)
                 //            - trans((F*S)*N())*idt(u)
                 + penalisation_bc*trans(idt(u))*id(v)/hFace() );
#endif
#if 0
    double Newmark_gamma=0.5;
    double Newmark_beta=1/4.;

    AUTO( cst_acc_displ, 1./(Newmark_beta*std::pow(M_bdf_displ_struct->timeStep(),2)) );

    form2( M_Xh, M_Xh, J ) +=
      integrate( elements(mesh), _Q<2*Order>(),
                 cst_acc_displ*trans(idt(u))*id(v) );
#endif
    J->close();

    //  J->addMatrix(1.0, this->OpLin );


  }

} // Life


// instantiation
template class Life::StVenantKirchhoff<2,2,Life::Simplex>;
//template class Life::StVenantKirchhoff<3,2,Life::Simplex>;

